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Abstract. Assemblies of granular particles mechanically stable under their own weight contain arches. 
These are structural units identified as sets of mutually stable grains. It is generally assumed that these 
arches shield the weight above them and should bear most of the stress in the system. We test such 
hypothesis by studying the stress born by in-arch and out-of-arch grains. We show that, indeed, particles 
in arches withstand larger stresses. In particular, the isotropic stress tends to be larger for in-arch-grains 
whereas the anisotropic component is marginally distinguishable between the two types of particles. The 
contact force distributions demonstrate that an exponential tail (compatible with the maximization of 
entropy under no extra constraints) is followed only by the out-of-arch contacts. In-arch contacts seem to 
be compatible with a Gaussian distribution consistent with a recently introduced approach that takes into 
account constraints imposed by the local force balance on grains. 
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1 Introduction 

The study of mechanically stable granular beds has be- 
come a major point of interest in granular Physics. Studies 
range from analysis of force network properties [IJ[21[21[1] to 
structural characterization [SlIHlEj to thermodynamic and 
statistical descriptions p'JTfll 1 Oil 1 T1fP2] . A recurrent ques- 
tion in the subject is related to the existence or not of a 
relation between force chains and arches [T51[T^] . 

Arches (or bridges) are multiparticle structures where 
all grains are mutually stable [13 , 15, 16, 17], i.e., fixing the 
positions of all other particles in the assembly, the removal 
of any particle in the arch leads to the destabilization of 
the other particles in it. For an arch to be formed, it is nec- 
essary (although not sufficient) that two or more falling 
particles be in contact at the time they reach mechanical 
equilibrium in order to create mutually stabilizing struc- 
tures [18] . Like arches in architecture, granular arches are 
assumed to sustain the weight of the material above. 

Highly stressed grains in static deposits are generally 
found to form linear structures: the so-called force chains. 
Notice that, the term "arch" is sometimes used [19,20,21 
to refer to these force chains and not to the mutually sta- 
bilizing structures defined above. Likewise, the term "dy- 
namic arch" has been used to refer to ephemeral structures 
that choke the flow of grains [35]. We have to distinguish 
between these usages and the classical meaning we adhere 
to in this work: an arch is a mechanically stable structure 
of mutually stabilizing bodies. 

Force chains are a clear spatial heterogeneity of the 
contact force network. The distribution of contact forces 



(both normal and tangential) shows no bimodal character. 
However, the spatial distribution of large and small forces 
is heterogeneous with large forces developing a somewhat 
open stringy network that encloses regions of grains that 
sustain little weight [23] , 

To what extent the bimodal spatial distribution of 
forces is related to the mutually stable structures (arches) 
has not been assessed so far. A correlation as been pointed 
out (13,15,24 between the distribution of horizontal span 
of arches in a granular pile and the distribution of normal 
forces at the grain contacts. Therefore, it is assumed that 
a strong correlation has to be present between arches and 
highly stressed grains in a granular deposit. In this pa- 
per, we assess this general belief in the frame of a simula- 
tion of granular packings prepared by tapping. The results 
provide additional information on the validity of basic as- 
sumptions made in the statistical description of granular 
packings. 



2 Simulation model 

To simulate packings of gains we have followed the stan- 
dard techniques on discrete element methods (see for ex- 
ample Refs. [25,26, 27J). We used a velocity Verlet algo- 
rithm to integrate the Newton equations for N monosized 
disks (diameter d) in a rectangular box of width L. We 
studied two system sizes: (i) N = 512 with L = 12.39c? and 
(ii) N = 2048 with L = 24.78d. The non-commensurate 
box is chosen to prevent crystallization to some extent. 
The larger system is roughly twice as tall as the smaller 
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system. However, this depends on the actual packing frac- 
tion obtained for a given preparation protocol. 

The disk-disk and disk-wall contact interaction com- 
prises a linear spring-dashpot in the normal direction 

F n = fc n C " 7n<, (1) 

and a tangential friction force 

F t = -min( M |F n |,|F s |)-sign(C) (2) 

that implements the Coulomb criterion to switch between 
dynamic and static friction [T%ll28| . 

In Eqs. Q-Q, £ = d — |ry| is the particle-particle 
overlap, represents the center-to-center vector between 
particles i and j, vfj is the relative normal velocity, F s = 

—k s (—^ s vl j is the static friction force, £ (i) = (t') dt' 

is the relative shear displacement, vj j = Yij-s+\d (uji + uij) 
is the relative tangential velocity, and s is a unit vector 
normal to ry. The shear displacement £ is calculated by 
integrating v\ ^ from the beginning of the contact (i.e., 
t = to). The disk-wall interaction is calculated consider- 
ing the wall as an infinite radius, infinite mass disk. Other 
than these, the interaction parameters are the same as for 
the disk-disk interaction. 

In these simulation we used the following set of pa- 
rameters: dynamic friction coefficient /i = 0.5, normal 
spring stiffness k n — 10 5 (mg/d), normal viscous damping 
7„ = 300(my/g/d), tangential spring stiffness k s = |fc n , 
and tangential viscous damping j s = 200 (to \/g~Jd). The 
integration time step is 5 — 10~ 4 Wd/g. Units are reduced 
with the diameter of the disks, d, the disk mass, to, and 
the acceleration of gravity, g. 

In order to investigate reproducible states, we imple- 
ment a tapping protocol. The system is initially deposited 
from a dilute configuration in which particles have no con- 
tacts nor overlaps. After the grains reach mechanical equi- 
librium, the system is tapped with a given amplitude and 
left to come back to mechanical equilibrium. After many 
taps of given amplitude, the system reaches a steady state 
where the properties of the static configurations gener- 
ated have well defined mean values and fluctuations. The 
steady state properties are independent of the initial con- 
ditions and are reproducible 12,29,31 . We generate a 
number of static packings after the steady state is reached 
to average quantities. Different steady states are generated 
by changing the tap amplitude. 

Tapping is simulated by moving the confining box in 
the vertical direction following a half sine wave trajec- 
tory of frequency v = n /2(g/d) 1/2 . The intensity of the 
excitation is controlled through the amplitude, A, of the 
sinusoidal trajectory; and it is characterized by the pa- 
rameter r = A{2irv) 2 j g . A new tap is applied only after 
the system has come to mechanical equilibrium, which is 
defined via the stability of each particle-particle contact 
[T5] , Averages were taken over 500 taps (configurations) 
in the steady state corresponding to each value of r after 
the 500 initial taps were discarded to avoid any transient. 



3 Identification of arches 

Details on the algorithms used to identify arches can be 
found in previous works 15,16,18 . Briefly, we need first 
to identify the supporting grains of each particle in the 
packing. In 2D, there are two disks that support any given 
grain. Two grains in contact with a given particle are able 
to provide support if the segment defined by the contact 
points lies below the center of mass of the particle. Some 
of these supporting contacts may be provided by the walls 
of the container. Then, we find all mutually stable parti- 
cles. Two grains A and B are mutually stable if A supports 
B and B supports A. Arches are defined as sets of parti- 
cles connected through these mutually stabilizing contacts 
(MSC). 

The fact that the supporting particles of each grain 
have to be known implies that contacts, and the chrono- 
logical order in which they form, have to be clearly de- 
fined in the model. Figure [I] shows some examples of the 
packings generated where arches are indicated by joining 
segments. In Fig. fljd), a close view of an arch detected in 
a given packing (formed by particles 1 to 5) is displayed. 
Particle 4 is an example of a grain whose pair of stabi- 
lizing particles is ambiguous form the limited information 
provided by the snapshot; discrimination requires chrono- 



(a) (b) (c) 




Fig. 1. (Color online). Sample images of the simulated gran- 
ular columns (N = 512) for different F: (a) r = 2.19, (b) 
r = 4.93 and (c) F = 15.39. The color code indicates the 
trace, Tr(cr), of the stress tensor for each particle in units of 
mg/d. The joining segments indicate the arches detected in the 
system, (d) A closeup on a 5-particle arch. See main text for a 
description on the supporting contacts of particle 4. 
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logical information. The pairs 3-5, 5-6 and 6-7 comply with 
the condition that the center of mass of grain 4 is above 
the segment that joins the corresponding contact points. 
However, the contacts with grains 3 and 5 where formed 
first and for that reason these particles are considered to 
be the supporting pair of grain 4. To identify the contacts 
that support each particle we use an algorithm that has 
been previously designed to work with molecular dynamic 
type simulations [18] . 



4 Single particle stress tensor 

We measure the stress tensor crj for grain i as [5]: 



aB 

r ~ 



7T(d/2) : 



.7=1 



b p - 



(3) 



where, /y is the force exerted by grain j on grain i and 
bij is the branch vector that goes from the center of grain 
i to the contact point with grain j. The sum runs over the 
N c particles in contact with particle i. 

The pressure (or isotropic stress) is given by the trace, 
Tr(cr), of a whereas the anisotropic component is charac- 
terized by the deviatoric stress s 




Fig. 2. (a) The mean packing fraction, <f>, as a function of 
tap intensity F. (b) Fraction of in-arch grains as a function of 
r. (c) The arch size distribution n(s) for F = 2.19. Results 
obtained for the two system sizes: N = 512 (solid triangles) 
and N = 2048 (open circles). 
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We use Dev(<r) = <t zz — cr xx to characterize the anisotropic 
component. In average, our packings under gravity present 
<j xx < a zz (i.e., the vertical component is higher than the 
horizontal component). 

The principal directions of the stress vary form config- 
uration to configuration during tapping. However, these 
fluctuations are very small since the shear component <r xz 
is less than 1% of Tr(er) in all our packings. 



5 General properties of the deposits 

In Fig. [2][a), we report the mean packing fraction, 0, as a 
function of the tap intensity, r. The packing fraction was 
measured in a slab of the bed that covers 50% of the height 
of the column and is vertically centered with the center of 
mass of the system. The values of </> for the smaller system 
is affected by the presence of the lateral walls, which tend 
to reduce the packing fraction. cj> presents a minimum at 
intermediate r as previously observed in various models 
[28, 3D] and experiments [T2]l31j . This minimum in <f> is 
related with the existence of a maximum in the number of 
grains involved in arches [see Fig. (2fb)] . A description of 
the mechanisms that lead to the existence of a maximum 
in the number of grains involved in arches can be found 
in Ref. [28] . 

In spite of the system being monodisperse, the pack- 
ings obtained present only partial crystallization. This is 
due to the non-commensurate simulation box. Even if very 
ordered packings were obtained, the contact forces (the 



main focus of this paper) have been found to display simi- 
lar statistics to the one shown by disordered packings [32 . 
We have also assessed the structural anisotropy through 
the fabric tensor. We found that all our packings present a 
deviatoric fabric of less than 5% of the fabric trace. There- 
fore, the structural anisotropy is rather small. 

The distribution, n(s), of the sizes of the arches found 
in the packings is shown in Fig. [2jc) . Here, n(s) is the 
fraction of arches consisting in s grains, with n(s = 1) the 
fraction of grains that do not belong to any arch. As we 
can see, n(s) is not affected by the system boundaries and 
arches of more than 10 disks have not been detected even 
in the 24-disk-wide system (i.e., N = 2048). 

Figure [T] shows some examples of the distribution of 
pressures and arches inside a granular pile. As it is to 
be expected, particles are subjected to higher pressures, 
in average, in the lower part of the pile as compared with 
the upper layers. The system does not display force chains 
that span the system from top to bottom as is commonly 
seen in many experiments and simulations. This is due 
to the fact that the system is in mechanical equilibrium 
under its own weight; no external compression is applied 
to the sample in any direction. 

In Fig. [3] we show the normal component of the contact 
forces for three different tap intensities (the lower and the 
higher r studied, and the value r min that leads to the 
minimum <p for the given system size). The mean normal 
contact force (F n ) increases rather linearly with the depth 
into the packing and is little dependent on the system size 
for any given depth [see Fig. |3][a)]. Only the system with 
2048 grains display a hint of Janssen saturation in the 
deeper layers. Small differences in (F n ) can be observed 
between packings obtained with different r. In particular, 
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for the lowest tap intensity considered, (F n ) is smaller at 
all depths. 

Figure |3jb) presents the normal contact force distri- 
bution for a depth of 35<i (this corresponds to the lower 
part of the smaller system and to the middle section of the 
larger system). All grain-grain contacts that lay within a 
slab 10<i-wide centered at a depth 35c? are considered. Tak- 
ing narrower slabs leads to similar results. As we can see, 
the PDF of F n coincides for both system sizes. We have 
seen that the tangential contact forces also show consis- 
tent results when systems of different sizes are compared 
by looking into slabs at the same depth. 

There exist a current debate on whether the tail of 
these distributions are or not exponentials [10,33 . Expo- 
nential tails for forces above the mean contact force have 
been reported by a number of authors considering granular 
packs subjected to external compression [2]J3"ll3"4] or stable 
under their own weight [IS]. As we can see in Fig. |3][b), 
for low r, we observe a clear exponential tail. However, a 
faster than exponential decay seem to be followed by the 
rest of the packings. Tighe et al. [TU] have argued that 
some reported exponential tails are perhaps Gaussians. 
The behavior for very small forces [see inset to Fig. |3jb)] 
resembles the weak divergence found for packings without 
external compression when bulk contacts (as opposed to 
contacts made between the grains and the container) are 
considered |3"4"1. 




V<F„> 

Fig. 3. (Color online), (a) Mean normal contact force (F n ) as a 
function of the depth into the granular column, (b) PDF of the 
normal contact force. We consider two system sizes: N = 512 
(solid symbols) and N — 2048 (open symbols). Results for 
three different tap intensities are reported (see legend). The 
intermediate value corresponds to the value of r that yields 
the minimum cj> for the given system size (i.e., r — 4.93 for 
N = 512 and T = 6.59 for N = 2048). The inset in part (b) is 
a close up for forces below the mean. 



In order to compare results from different system sizes, 
the remaining of the paper, unless otherwise specified, will 
refer to measurements made in a slab lOd-wide centered 
at a depth 35d. 



6 Contact forces and arches 

It can be observed from Fig. [T] that, at any depth into the 
pile, grains can present high and low stress irrespective of 
whether they belong to an arch or not. Also, force chains 
do not coincide with arches although arches form part of 
portions of these chains. For a more quantitative analy- 
sis we plot in Fig. [4] the mean value of the contact forces 
(normal and tangential to the contact in a slab at 35d of 
depth) as a function of r. MSC (mutually stabilizing con- 
tacts) and non-MSC have been separated in the analysis. 
Although some small differences are observed between the 
results for the two system sizes studied, the general trends 
are quite similar. It is clear that MSC (i.e., contacts within 
arches) have, in average, larger (roughly 50%) normal and 
tangential forces. This supports the idea that arches bear 
most of the stress in the system and that force chains and 
arches must be correlated. 

Figure[5ja) shows that the distribution of contact forces 
for MSC and non-MSC are clearly distinguishable for the 
normal component. Although, we present the distribution 
obtained for F m i n , most packings display the same general 
features (some exceptions regarding packings prepared at 
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Fig. 4. (Color online). The mean value of the contact force for 
mutually stabilizing contacts (blue) and non-mutually stabiliz- 
ing contacts (red), (a) Normal contact forces, and (b) tangen- 
tial contact forces. Results obtained for the two system sizes: 
N = 512 (solid triangles) and N = 2048 (open circles). 
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Fig. 5. (Color online). The PDF of contact forces for MSC 
(blue) and non-MSC (red) for _T = r min . (a) Normal contact 
forces, and (b) tangential contact forces. Results obtained for 
the two system sizes: N = 512 (solid triangles) and N — 2048 
(open circles). 



low r are discussed below). The mild divergence for very 
small forces is still present for both distributions. Despite 
the difference, there is not a clear separation of the two 
populations of contacts. The bimodal character observed 
in the spatial distribution of contacts seems to be poorly 
correlated with MSC. 

As we can see from Fig. [5|a), the PDF for non-MSC 
present a clear exponential tail, whereas MSC present a 
faster decay. The non-MSC PDF corresponds to an expo- 
nential decay even for forces below the mean. The expo- 
nential PDF has a well established statistical explanation. 
If the mean force is set and all contact force states are 
equally probable, an exponential distribution of contact 
forces maximizes the entropy (defined as the logarithm of 
the number of contact force states) [361135] . 

MSC seem to have a distribution compatible with a 
Gaussian tail, or at least a faster than exponential tail. 
It seems that the deviation form an exponential in the 
full PDFs reported in the literature [and in Fig. [3^b)] 
seem to be due to the presence of MSC (and therefore 
the presence of arches) . The immediate conclusion is that 
the presence of arches prevents us from making some of 
the basic assumptions on the contact forces to render a 
simplified statistical analysis. In particular, arches intro- 
duce force balance constraints that need to be accounted 
for. Tighe et al. [TU] have shown that force balance con- 
straints (the conservation of the total area of the Maxwell 
reciprocal tiling) can be introduced in a force ensemble. 
These have led to Gaussian contact force distributions. 
Notice however that this theoretical approach yields the 



same Gaussian distribution irrespective of the existence of 
arches in the packing. 

Figure [fjjb) shows that the tangential components of 
the contact forces have a much subtle difference between 
the distribution for MSC and non-MSC. Again, MSC present 
a somewhat faster-than-exponential tail in contrast with 
the non-MSC. 

We now focus in the results for the smallest tap in- 
tensity reported. As we mentioned, Fig. [3j^b) shows that 
for r = 2.19 the PDF for normal contacts presents a 
clear exponential tail, in contrast with the packings gen- 
erated with stronger taps. Separating MSC and non-MSC 
in the analysis leads to two exponential tails (presenting 
slightly different slopes). We speculated that there could 
be fewer MSC in these packings than in packings obtained 
with stronger taps. However, these packings present simi- 
lar number of MSC as compared with packings that show a 
faster-than-exponential tail. The main difference we have 
been able to find is that these packings have, in compar- 
ison, fewer arches composed of three or more grains. It 
seems that arches composed of three or more particles 
are the responsible for introducing strong force balance 
constraints that render the PDF non-exponential. Some 
reports of pure exponential decays can be found in previ- 
ous studies. Blair et al. mentioned that a pure exponential 
was found in some cases depending on the history of the 
packing [35]. Makse et al. found pure exponentials too 
in simulations of isotropically compressed grains, which 
may develop structures without arches [37]. We believe 
the preparation history of these packings may have lead 
to a small presence of arches composed of three or more 
grains. 

7 Stress tensor and arches 

In Fig. [6] we show the results of an analysis similar to the 
previous section but now the stress tensor on each parti- 
cle, as defined in Eq. ([3]), is considered. The stress tensor 
accounts for both MSC and non-MSC on each grain. We 
separate in-arch grains from out-of-arch grains in the anal- 
ysis. In-arch grains support, in average, isotropic pressures 
[see Tr(er) in Fig. (6fa)] about 20% higher than out-of- 
arch grains. In contrast, the anisotropic component of the 
stress, Dev(cr), seems to be rather similar for both types of 
grains. This implies that the actual difference between the 
stress tensor of in-arch and out-of-arch grains corresponds 
to the addition of a constant to the diagonal components 
(in contrast to an increse given by a multiplicative con- 
stant). We have seen that the shear stress a xz is the same 
for both types of grains. 

Figure]^ a) shows that the distribution of the isotropic 
stress is markedly different for in-arch and out-of-arch 
grains. In-arch grains present a clear maximum in the 
PDF of Tr(cr) at around the mean. Although there is not a 
strong separation, the maximum in the PDFs suggest that 
the well known bimodal character of the force network is 
driven by the presence of arches to some extent. 

The distribution of the stress deviator is presented in 
Fig. [7][b). As we can see, the PDFs for in-arch and out-of- 
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Fig. 6. (Color online). Mean stress tensor for in-arch (blue) 
and out-of-arch (red) grains, (a) The trace Tr(cr) of the stress, 
and (b) the deviator Dev(rj). Results obtained for the two sys- 
tem sizes: N = 512 (solid triangles) and N — 2048 (open 
circles) . 



arch grains are almost identical. The negative values are 
due to the fact that some grains have a xx > a yy . However, 
in average a xx < a yy and the mean deviator as defined 
above is always positive in our packings. 



8 Conclusions 

We have shown that MSC, which define arches, present 
higher normal and tangential components of the contact 
forces as compared with non-MSC. Grains belonging to 
arches are generally subjected to larger isotropic stresses 
but similar anisotropic stress. Therefore, particles in arches 
are, to some extent, different from particles that do not 
form arches when their contact forces are considered. This 
is in line with the common assumption that arches carry 
most of the weight in a granular deposit. 

The PDF of normal contact forces show that non-MSC 
follow an exponential decay whereas the MSC present a 
faster-than-exponential fall. This has strong implications 
for the basic statistical models of force distribution. In 
particular, it seems that MSC are the main cause for the 
constraints in force balance not considered in simplistic 
approaches. These constraints lead to the deviation of the 
overall-contacts PDF from the expected exponential. In- 
deed, packings containing a low number of large arches 
(arches of three or more grains) seem to fit better the 
exponential law. 

The bimodal spatial distribution of stresses seems to 
be related to some extent with the presence of arches. 
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Fig. 7. (Color online). The PDF of the stress tensor for in- 
arch (blue) and out-of-arch (red) for r = r min . (a) The trace, 
Tr(er), of the stress tensor, and (b) the deviator Dev(tr). Results 
obtained for the two system sizes: N = 512 (solid triangles) and 
N = 2048 (open circles). 



Particles in arches present a clear maximum around the 
mean stress in the PDF of isotropic stress. 

It is worth mentioning that despite the correlations 
found between arches and force chains, there is not a one- 
to-one correspondence. Arches that sustain little weight 
can always exist in the structure since they are shielded 
by other arches above. This leads to the preponderance of 
very small forces in the distributions for MSC. Also, force 
chains can develop without the need of arches. A deposit 
carefully built by sequential deposition of grains contains 
no arches in the structure, yet it will present force chains. 
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